
#######################################################
#######################################################

# Replication code for "Oil Income and the Personalization of Autocratic Politics"
# Author: Matthew Fails
# Email: 	fails@oakland.edu
#PSRM manusctipt #PSRM-RN-2018-0133.R1

## R script for Appendix Table 3A, Robustness to incorporation of uncertainty in latent dependent variable
## This script contains the code to reproduce Table 3A Models 7 and 8. Each reads in a separate csv datafile, processed outside of R to match specifications in main results
## Script follows Crabtree and Fariss (2015) and Schnakenberg and Fariss (2014)

#######################################################
#######################################################

## Table 3A, model 7 - post Cold War sample, with additional control variables

##################

# setwd("C:/data")
# getwd()

##################

library(plm)
library(foreign)
library(sandwich)
library(lmtest)

#Read in the dataset

data <- read.csv("AppModel7.csv", header=TRUE)
head(data)

#Create a list to save datasets
newdata <- list()

################

datnew <- data

################

# Take K draws from the posterior distribution and make K new datasets in the list object just created

set.seed(123456789)

for(i in 1:1000){
  newdata[[i]] <- datnew
  newdata[[i]]$pers_2pl <- rnorm(cbind(rep(1,nrow(data))), mean=datnew$pers_2pl, sd=datnew$pers_se_2pl)

}
#######################

FML <- "pers_2pl ~ lag2oil + lag2inc + lagpop + lagnew + lagdur + lagfdi + lagtax + yrdum1 + yrdum2 + yrdum3 + yrdum4 + yrdum5 + yrdum6 + yrdum7 + yrdum8 + yrdum9 + yrdum10 + yrdum11 + yrdum12 + yrdum13 + yrdum14 + yrdum15 + yrdum16 + yrdum17 + yrdum18 + caseid1 + caseid2 +	caseid3	+	caseid4	+	caseid5	+	caseid6	+	caseid7	+	caseid8	+	caseid9	+	caseid10	+	caseid11	+	caseid12	+	caseid13	+	caseid14	+	caseid15	+	caseid16	+	caseid17	+	caseid18	+	caseid19	+	caseid20	+	caseid21	+	caseid22	+	caseid23	+	caseid24	+	caseid25	+	caseid26	+	caseid27	+	caseid28	+	caseid29	+	caseid30	+	caseid31	+	caseid32	+	caseid33	+	caseid34	+	caseid35	+	caseid36	+	caseid37	+	caseid38	+	caseid39	+	caseid40	+	caseid41	+	caseid42	+	caseid43	+	caseid44	+	caseid45	+	caseid46	+	caseid47	+	caseid48	+	caseid49	+	caseid50	+	caseid51	+	caseid52	+	caseid53	+	caseid54	+	caseid55	+	caseid56	+	caseid57	+	caseid58	+ caseid59	+	caseid60	+	caseid61	+	caseid62	+	caseid63	+	caseid64 + caseid65 + caseid66	+	caseid67	+	caseid68	+	caseid69	+	caseid70	+	caseid71	+	caseid72	+	caseid73	+	caseid74	+	caseid75	+	caseid76	+	caseid77	+	caseid78	+	caseid79	+	caseid80	+	caseid81	+	caseid82"

#########################
# Create function. The code below runs a linear model K times and combines all the estimates using the formula for combining multiply imputed datasets from Rubin.
milm <- function(fml, midata){
  xx <- terms(as.formula(fml))
  lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcovHC(tmp, type = "HC")))
    vcovs[[i]] <- vcovHC(tmp, type = "HC")
  }
  par.est <- apply(lms, 1, mean)
  se.within <- apply(ses, 1, mean)
  se.between <- apply(lms, 1, var)
  se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
  list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)   
}

# Call function
model <- milm(fml=FML, midata=newdata)

# Inspect returned object
attributes(model)

# Organize results
tstats <- data.frame(model$terms,model$beta,model$SE,model$beta/model$SE)
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat")

# View results
tstats


#########################################################
#########################################################

## Table 3A, model 8 - post Cold War sample, with additional control variables, exclude top producers

##################

# setwd("C:/data")
# getwd()

##################

library(plm)
library(foreign)
library(sandwich)
library(lmtest)

#Read in the dataset

data <- read.csv("AppModel8.csv", header=TRUE)
head(data)

#Create a list to save datasets
newdata <- list()

################

datnew <- data

################

# Take K draws from the posterior distribution and make K new datasets in the list object just created

set.seed(123456789)
	
for(i in 1:1000){
  newdata[[i]] <- datnew
  newdata[[i]]$pers_2pl <- rnorm(cbind(rep(1,nrow(data))), mean=datnew$pers_2pl, sd=datnew$pers_se_2pl)

}
#######################

FML <- "pers_2pl ~ lag2oil + lag2inc + lagpop + lagnew + lagdur + lagfdi + lagtax + yrdum1 + yrdum2 + yrdum3 + yrdum4 + yrdum5 + yrdum6 + yrdum7 + yrdum8 + yrdum9 + yrdum10 + yrdum11 + yrdum12 + yrdum13 + yrdum14 + yrdum15 + yrdum16 + yrdum17 + yrdum18 + caseid1 + caseid2 +	caseid3	+	caseid4	+	caseid5	+	caseid6	+	caseid7	+	caseid8	+	caseid9	+	caseid10	+	caseid11	+	caseid12	+	caseid13	+	caseid14	+	caseid15	+	caseid16	+	caseid17	+	caseid18	+	caseid19	+	caseid20	+	caseid21	+	caseid22	+	caseid23	+	caseid24	+	caseid25	+	caseid26	+	caseid27	+	caseid28	+	caseid29	+	caseid30	+	caseid31	+	caseid32	+	caseid33	+	caseid34	+	caseid35	+	caseid36	+	caseid37	+	caseid38	+	caseid39	+	caseid40	+	caseid41	+	caseid42	+	caseid43	+	caseid44	+	caseid45	+	caseid46	+	caseid47	+	caseid48	+	caseid49	+	caseid50	+	caseid51	+	caseid52	+	caseid53	+	caseid54	+	caseid55	+	caseid56	+	caseid57	+	caseid58	+ caseid59	+	caseid60	+	caseid61	+	caseid62	+	caseid63	+	caseid64 + caseid65 + caseid66	+	caseid67	+	caseid68	+	caseid69	+	caseid70	+	caseid71	+	caseid72 + caseid73 + caseid74 + caseid75 + caseid76 + caseid77 + caseid78 + caseid79 + caseid80 + caseid81 + caseid82"

#########################

milm <- function(fml, midata){
  xx <- terms(as.formula(fml))
  lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcovHC(tmp, type = "HC")))
    vcovs[[i]] <- vcovHC(tmp, type = "HC")
  }
  par.est <- apply(lms, 1, mean)
  se.within <- apply(ses, 1, mean)
  se.between <- apply(lms, 1, var)
  se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
  list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)   
}

# Call function
model <- milm(fml=FML, midata=newdata)

# Inspect returned object
attributes(model)

# Organize results
tstats <- data.frame(model$terms,model$beta,model$SE,model$beta/model$SE)
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat")

# View results
tstats
